//============================================================================
// Name        : lagrange-spline.cpp
// Author      : Tomasz Lipert
// Version     :
// Copyright   : 
// Description : Obliczanie wartości i współczynników naturalnej funkcji sklejanej
//				 stopnia trzeciego oraz wartości i współczynników wielomianu interpolacyjnego
//				 Lagrange'a - porównanie.
//============================================================================

#include <iostream>
using namespace std;

int n = 6;
long double jaki_x;
int st;

//testowe dane	cout <<
long double x[] = { 17, 20, 23, 24, 25, 27, 27.7 };
long double y_in[] = { 4.5, 7, 6.1, 5.6, 5.8, 5.2, 4.1 };
long double a[10][10];

//long double x[] = { 0, 1 };
//long double y_in[] = { 0, 1 };


void wczytaj_dane() {
	cout << "Podaj liczbę węzłów: ";
	cin >> n;

	for (int i = 0; i < n; i++) {
		cout << "Podaj x[" << i << "]: ";
		cin >> x[i];
		cout << "Podaj y[" << i << "]: ";
		cin >> y_in[i];

	}

}

//funkcja wrażliwa na kolejność węzłów
long double lagrangeInterpolation(long double xs[], long double ys[],
		long double x) {
	long double t;
	long double y = 0.0;

	for (int k = 0; k < n; k++) {
		t = 1.0;
		for (int j = 0; j < n; j++) {
			if (j != k) {
				t = t * ((x - xs[j]) / (xs[k] - xs[j]));
			}
		}
		y += t * ys[k];
	}
	return y;
}

long double naturalsplinevalue(int n, long double x[], long double f[],
		long double xx) {

	int i, k, st;
	long double u, y, z;
	bool found;
	long double a[4] = { 0, 0, 0, 0 };
	long double d[7] = { 0 }, b[7] = { 0 }, c[7] = { 0 };

	if (n < 1)
		st = 1;
	else if ((xx < x[0]) || (xx > x[n]))
		st = 3;
	else {
		st = 0;
		i = -1;
		do {
			i = i + 1;
			for (k = (i + 1); k <= n; k++) {
				if (x[i] == x[k]) {
					st = 2;
				}
			}
		} while (!(i == n - 1) && !(st == 2));
	}
	if (st == 0) {
		d[0] = 0;
		d[n] = 0;
		if (n > 1) {
			z = x[2];
			y = z - x[1];
			z = z - x[0];
			u = f[1];
			b[1] = y / z;
			d[1] = 6 * ((f[2] - u) / y - (u - f[0]) / (x[1] - x[0])) / z;
			z = x[n - 2];
			y = x[n - 1] - z;
			z = x[n] - z;
			u = f[n - 1];
			c[n - 1] = y / z;
			d[n - 1] = 6
					* ((f[n] - u) / (x[n] - x[n - 1]) - (u - f[n - 2]) / y) / z;
			int indeks = 2;
			for (indeks = 2; indeks <= (n - 2); indeks++) {
				z = x[indeks];
				y = x[indeks + 1] - z;
				z = z - x[indeks - 1];
				u = f[indeks];
				b[indeks] = y / (y + z);
				c[indeks] = 1 - b[indeks];
				d[indeks] = 6 * ((f[indeks + 1] - u) / y - (u - f[indeks - 1])
						/ z) / (y + z);
			}
			u = 2;
			i = 0;
			y = d[1] / u;
			d[1] = y;
			if (n > 2) {
				do {
					i = i + 1;
					z = b[i] / u;
					b[i] = z;
					u = 2 - z * c[i + 1];
					y = (d[i + 1] - y * c[i + 1]) / u;
					d[i + 1] = y;
				} while (!(i == n - 2));
			}
		}
		u = d[n - 1];
		for (i = n - 2; i >= 1; i--) {

			u = d[i] - u * b[i];
			d[i] = u;
		}
		found = 0;
		i = -1;
		do {
			i = i + 1;
			if ((xx >= x[i]) && (xx <= x[i + 1]))
				found = 1;
		} while (!(found));
		y = x[i + 1] - x[i];
		z = d[i + 1];
		u = d[i];
		a[0] = f[i];
		a[1] = (f[i + 1] - f[i]) / y - (2 * u + z) * y / 6;
		a[2] = u / 2;
		a[3] = (z - u) / (6 * y);
		y = a[3];
		z = xx - x[i];
		for (i = 2; i >= 0; i--) {
			y = y * z + a[i];
		}
	}
	if (st == 0) {
		return y;
	} else
		cout << "Wystąpił wyjątek. ST jest równe " << st << "     ";
}

void naturalsplinecoeffns(int n, long double x[], long double f[]) {
	int i, k;
	long double u, v, y, z, xi;
	long double d[7] = { 0 }, b[7] = { 0 }, c[7] = { 0 };

	if (n < 1)
		st = 1;
	else {
		st = 0;
		i = -1;
		do {
			i = i + 1;
			for (k = i + 1; k <= n; k++) {
				if (x[i] == x[k])
					st = 2;
			}
		} while (!(i == n - 1) && !(st == 2));
	}
	if (st == 0) {
		d[0] = 0;
		d[n] = 0;
		if (n > 1) {
			z = x[2];
			y = z - x[1];
			z = z - x[0];
			u = f[1];
			b[1] = y / z;
			d[1] = 6 * ((f[2] - u) / y - (u - f[0]) / (x[1] - x[0])) / z;
			z = x[n - 2];
			y = x[n - 1] - z;
			z = x[n] - z;
			u = f[n - 1];
			c[n - 1] = y / z;
			d[n - 1] = 6
					* ((f[n] - u) / (x[n] - x[n - 1]) - (u - f[n - 2]) / y) / z;
			for (i = 2; i <= n - 2; i++) {
				z = x[i];
				y = x[i + 1] - z;
				z = z - x[i - 1];
				u = f[i];
				b[i] = y / (y + z);
				c[i] = 1 - b[i];
				d[i] = 6 * ((f[i + 1] - u) / y - (u - f[i - 1]) / z) / (y + z);
			}
			u = 2;
			i = 0;
			y = d[1] / u;
			d[1] = y;
			if (n > 2) {
				do {
					i = i + 1;
					z = b[i] / u;
					b[i] = z;
					u = 2 - z * c[i + 1];
					y = (d[i + 1] - y * c[i + 1]) / u;
					d[i + 1] = y;

				} while (!(i == n - 2));
			}
			u = d[n - 1];
			for (i = n - 2; i >= 1; i--) {

				u = d[i] - u * b[i];
				d[i] = u;
			}
			for (i = 0; i <= n - 1; i++) {
				u = f[i];
				xi = x[i];
				z = x[i + 1] - xi;
				y = d[i];
				v = (f[i + 1] - u) / z - (2 * y + d[i + 1]) * z / 6;
				z = (d[i + 1] - y) / (6 * z);
				y = y / 2;
				a[0][i] = ((-z * xi + y) * xi - v) * xi + u;
				u = 3 * z * xi;
				a[1][i] = (u - 2 * y) * xi + v;
				a[2][i] = y - u;
				a[3][i] = z;
			}
		}
	}
}

void licz_wyswietl_lagrange() {
	//	for (jaki_x = 0; jaki_x < 10; jaki_x = jaki_x + 0.1) {
	long double y2 = lagrangeInterpolation(x, y_in, jaki_x);
	//	cout << "x = " << jaki_x << " -> y = " << y2 << "\n";
	cout << jaki_x << "  " << y2 << "\n";
	//}
}

void licz_wyswietl_spline() {
	//for (jaki_x = 0; jaki_x < 10; jaki_x = jaki_x + 0.1) {
	jaki_x = 23.5;
	long double y2 = naturalsplinevalue(n, x, y_in, jaki_x);
	//cout << "x = " << jaki_x << " -> y = " << y2 << "\n";
	cout << jaki_x << "  " << y2 << "\n";
	//}
}

void licz_wyswietl_spline_wsp() {
	naturalsplinecoeffns(n, x, y_in);
	for (int i = 0; i < 10; i++) {
		for (int j = 0; j < 10; j++) {
			cout << a[j][i] << "    ";
		}
		cout << "\n";
	}
}

int main() {
	//wczytaj_dane();
	//cout << "Lagrange\n";
	//licz_wyswietl_lagrange();
	cout << "Splajn\n";
	cout.precision(5);
	//licz_wyswietl_spline();
	licz_wyswietl_spline_wsp();
	return 0;
}
